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This paper presents a rewriting-logic -based modeling and analysis technique for physical systems, 
with focus on thermal systems. The contributions of this paper can be summarized as follows: (i) 
providing a framework for modeling and executing physical systems, where both the physical compo- 
nents and their physical interactions are treated as first-class citizens; (ii) showing how heat transfer 
problems in thermal systems can be modeled in Real-Time Maude; (iii) giving the implementation 
in Real-Time Maude of a basic numerical technique for executing continuous behaviors in object- 
oriented hybrid systems; and (iv) illustrating these techniques with a set of incremental case studies 
using realistic physical parameters, with examples of simulation and model checking analyses. 

1 Introduction 

The rewriting-logic-based Real-Time Maude tool lfl2l has proved to be very useful for formally modeling 
and analyzing a wide range of advanced real-time systems (see, e.g., J6j|9j[lll[l3j[l4l) that are beyond 
the scope of timed- automaton-based tools. One important question to investigate is to what degree - and 
how - Real-Time Maude can be successfully applied to formally model and analyze hybrid systems. 

As part of this investigation, this paper presents a rewriting-logic-based modeling and analysis tech- 
nique for physical systems, which consist of a set of components that behave and interact according to the 
laws of physics. On the one hand, the physical quantities of the components continuously evolve in time. 
On the other hand, the components may also have discrete state parts, thus exhibiting hybrid behaviors. 

We propose a technique to generate executable models of physical systems in Real-Time Maude. 
Our modeling technique is based on an adaptation of the effort and flow approach described in ifTTl . 
explicitly modeling the transfer of power between the components. For the continuous behavior of 
physical systems, which is described by differential equations, the execution is based on the Euler method 
(see, e.g., Q) giving approximate solutions to ordinary differential equations. This numerical method 
operates with a discrete sampling strategy over a continuous time domain. 

This paper focuses on thermal systems, which are physical systems that deal with heat transfer 
and temperature change problems. Thermal systems appear in many important computer-controlled 
(or computer-controllable) applications, such as house heating systems, human body thermoregulatory 
systems, and nuclear power plants. We describe the use of our modeling and execution technique with a 
sequence of increasingly complex thermal systems. We first model simple thermal systems, such as "a 
cup of coffee in a room"; then we add a coffee heater to the system; and, finally, we add an automatic 
coffee heater as a controlling system that manages the coffee temperature. 

The main contributions of this paper can be summarized as follows: 

• We develop a framework for modeling general physical systems in Real-Time Maude. 

• This framework is used as the basis to define a technique to model heat transfer in thermal systems. 
We provide some basic classes, equations, and rules that can be used or extended to build models 
of thermal systems. 
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• The implementation of a basic numerical technique to execute continuous behavior of object- 
oriented real-time systems, which allows the analysis of their hybrid behavior. 

• Incremental case studies, using realistic physical parameters, with examples of simulation and 
reachability analysis. 

We give an overview of Real-Time Maude in Section|2| Section [3] gives a brief introduction to ther- 
mal systems. Section|4]presents a high-level overview of modeling physical and thermal systems with the 
effort and flow approach, and approximating their behavior. Section|5]describes our modeling and execu- 
tion framework for thermal systems using Real-Time Maude. Several case studies are given in Section[6] 
Section [TJdiscusses related work, and Section [8] mentions future work and summarizes the paper. 

2 Real-Time Maude 

In this section we briefly introduce Real-Time Maude lfT2l , a rewriting-logic-based tool supporting the 
formal modeling and analysis of real-time systems. A Real-Time Maude timed module specifies a real- 
time rewrite theory (L,E,IR,TR), where: 

• (L,E) is a membership equational logic JU theory with £ a signatureQand E a set of confluent and 
terminating conditional equations. specifies the system's state space as an algebraic data 
type, and must contain a specification of a sort Time modeling the (discrete or dense) time domain. 

• IR is a set of (possibly conditional) labeled instantaneous rewrite rules of the form 

rl [/] : t => t' 
crl [Z] : t => t 1 if cond 

specifying the system's instantaneous (i.e., zero-time) one-step transitions from an instance of t to 
the corresponding instance of t' , where I is a label. Conditional instantaneous rules apply only if 
their conditions hold. The rules are applied modulo the equations E^ 

• TR is a set of (possibly conditional) tick (rewrite) rules, written with syntax 

rl [/] : {t} => {/} in time T 
crl Ul : {?} => {?'} in time X if cond 

that model time elapse. {_} is a built-in constructor of sort GlobalSystem, and x is a term of sort 
Time that denotes the duration of the rewrite. 

The initial state must be a ground term of sort GlobalSystem and must be reducible to a term of the form 
{?} using the equations in the specifications. The form of the tick rule then ensures that time advances 
uniformly in all parts of the system. 

The Real-Time Maude syntax is fairly intuitive. For example, a function symbol / is declared with 
the syntax op / : s\ . . . s n -> s, where s\ ... s n are the sorts of its arguments, and s is its (value) sort. 
Equations are written with syntax eq t = t', and ceq t = t' if cond for conditional equations. The 
mathematical variables in such statements are declared with the keywords var and vars. We refer to El 
for more details on the syntax of Real-Time Maude. 

In object-oriented Real-Time Maude modules, a class declaration 

'i.e., L is a set of declarations of sorts, subsorts, and function symbols 

2 E is a union E 1 UA, where A is a set of equational axioms such as associativity, commutativity, and identity, so that 
deduction is performed modulo A. Operationally, a term is reduced to its E'-normal form modulo A before any rewrite rule is 
applied. 
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class C I att\ : s\, ... , att„ : s n . 

declares a class C with attributes att\ to att n of sorts s\ to s„. An object of class C in a given state 
is represented as a term < : C | atti : val\, ...,att n : va/ n > of sort Object, where O, of sort Oid, is 
the object's identifier, and where val\ to val„ are the current values of the attributes att\ to att n . In a 
concurrent object-oriented system, the state is a term of the sort Configuration. It has the structure of 
a multiset made up of objects and messages (that are terms of sort Msg). Multiset union for configurations 
is denoted by a juxtaposition operator (empty syntax) that is declared associative and commutative, so 
that rewriting is multiset rewriting supported directly in Real-Time Maude. 

The dynamic behavior of concurrent object systems is axiomatized by specifying its concurrent tran- 
sition patterns by rewrite rules. For example, the rule 

rl [1] : < : C I al : x, a2 : y, a3 : z > 

< 0' : C I al : w, a2 : 0, a3 : v > 
=> 

< : C I al : x + w, a2 : y, a3 : z > 

< 0' : C I al : w, a2 : x, a3 : v > 

defines a family of transitions where two objects of class C synchronize to update their attributes when 
the a2 attribute of one of the objects has value 0. The transitions have the effect of altering the attribute 
al of the object and the attribute a2 of the object ' . "Irrelevant" attributes (such as a3 and a2 of 0, 
and the right-hand side occurrence of al of ') need not be mentioned in a rule. 
A subclass inherits all the attributes and rules of its superclasses. 

A Real-Time Maude specification is executable under reasonable conditions, and the tool offers a 
variety of formal analysis methods. The rewrite command simulates one fair behavior of the system up 
to a certain duration. It is written with syntax (trew t in time <= % . ), where t is the initial state 
and T is a term of sort Time. The search command uses a breadth-first strategy to analyze all possible 
behaviors of the system, by checking whether a state matching a pattern can be reached from the initial 
state such that a given condition is satisfied. The timed search command, having the syntax 

(tsearch [1] t =>* pattern such that cond in time <= T . ) 

works similarly, but restricts the search to states reachable from the initial state within time t. 

Real-Time Maude also extends Maude's linear temporal logic model checker to check whether each 
behavior, possibly up to a certain time bound, satisfies a linear temporal logic formula. State propositions 
are terms of sort Prop, and their semantics should be given by (possibly conditional) equations of the 
form istatePattern} I = prop = b, for b a term of sort Bool, which defines the state proposition prop to 
hold in all states it} where it} \ = prop evaluates to true. A temporal logic formula is constructed by 
state propositions and temporal logic operators such as True, False, ~ (negation), A, \/, ~> (implica- 
tion), [] ("always"), <> ("eventually"), and U ("until"). The time-bounded model checking command 
has syntax 

(mc t |=t formula in time <= T . ) 

for initial state t and temporal logic formula/ormw/a . 

Finally, the find earliest command determines the shortest time needed to reach a desired state. 

For time-nondeterministic tick rules (i.e., tick rules in which the matching substitution does not 
uniquely determine the duration of the tick), the above model checking commands are applied according 
to the chosen time sampling strategy, so that only a subset of all possible behaviors is analyzed. 
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3 Thermal Systems 

We define a thermal system as a system whose components behave according to the laws of physics and 
that are related to heat transfer and temperature change. Heat is the form of energy that can be transferred 
from one system to another as a result of temperature difference. The transfer of energy as heat is always 
from the higher-temperature medium to the lower-temperature one. The flow of heat to a component 
effects the component by either changing its temperature or by changing the component's phase (such as 
ice melting to water). 

Temperature change. The equation representing the heat gained or lost by an object as it undergoes a 
temperature change AT is given by AQ = m-c- AT, where m is the mass of the object and c is the specific 
heat of the object at constant volume. The amount of heat transferred per time unit, the heat transfer 
rate, can be described as Q = m ■ c • t , where t represents the change of temperature per time unit. 
Heat can be transferred by three different mechanisms: conduction, convection, and radiation. 

• Conduction occurs when heat flows through stationary materials. The rate of heat conduction 
through a medium depends on the geometry of the medium, the thickness of the medium, the 
material of the medium, and temperature difference across the medium. The conduction rate is 
given by Q = ^ • (7\ — T2), where k is the thermal conductivity of the material, A is the area of 
the conduction, L is the thickness of the material through which the conduction occurs, and T\ , Ti 
are the temperatures of the two different media. 

• Convection occurs when a moving fluid transports heat from a hotter object to a colder object. The 
rate of heat convection is proportional to the temperature difference, and is expressed by Newton's 
law of cooling as Q = h ■ A ■ (T\ — T2), where h is the convection heat transfer coefficient, and A is 
the surface area through which convection heat transfer takes place. 

• Radiation occurs directly through the components' surface, without any transfer medium. The rate 
of heat radiation is given by Q = e ■ a ■ A ■ (T* — T 2 4 ), where £ is the emissivity of the surface, a 
is the Stefan-Boltzmann constant, and A is the surface area through which radiation heat transfer 
takes place. 

For a running example, consider a cup of hot coffee in a room. If the temperature of the coffee is 
higher than that of the room, heat will flow from the coffee to the room, until both of them reach the 
same temperature. In particular, the heat flows through the cup wall by conduction and radiation, and 
from the surface of the coffee (assuming there is no lid on the cup) to the room by convection. 

Phase transition. Another important aspect that must be taken into account when modeling thermal 
systems is the different phases of a substance. For example, water has three distinct phases: solid (ice), 
liquid, and gas (vapor). The latent heat is the amount of energy released or absorbed by a chemical 
substance during a phase transition. Note that during a phase transition the temperature does not change; 
the flow of energy goes in as the latent heat of the transition process. For example, to change its phase 
from solid to liquid, water goes through the melting process which consumes some energy. 

4 A Framework for Modeling Physical Systems 

This section presents a high-level overview of the general method we use to model physical systems 
and to approximate their behavior. This method adapts the effort-and-flow-variable approach common 
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Figure 1: Physical system components. 



in the field of physical systems modeling (see, e.g., 1171 ). In particular, Section 4.1 briefly describes 
the modeling of physical systems in this framework, and Section |4.2| shows how our method can be 
instantiated to model thermal systems. 



4.1 Modeling Physical Systems 

A well known modeling approach for physical systems is based on two kinds of variables in the model 
specification: effort and flow variables ifTTl . 

The components of a physical system can be thought of as energy manipulators which process the 
energy injected into the system. These energy manipulators interact through energy ports. For this 
purpose, effort and flow variables are used to represent the power being transmitted through energy 
ports. The flow variable is associated with the act of delivering energy, whereas the effort variable is 
associated with the act of measuring the flow of energy. For example, our coffee system could have 
two effort variables, one denoting the temperature of the coffee and one denoting the temperature of the 
room, and one flow variable, denoting the flow of heat from the coffee to the room per time unit. 

The approach using effort/flow variables is applicable to different areas of physical systems. In me- 
chanical translation systems, the pair of effort and flow variables are force and velocity; in mechanical 
rotation systems, torque and angular velocity; in electrical systems, voltage and current; in fluidic sys- 
tems, pressure and volume flow rate; in thermal systems, temperature and heat flow rate IfTTl . 

Our idea is to think of a physical system as the combination of two kinds of components: physical 
entities and physical interactions, as shown in Fig.[T] 

A physical entity (such as a cup of coffee) consists of a set of attributes and a continuous dynamics. 
We consider three kinds of attributes: continuous variables (denoting physical quantities, such as the 
temperature, that change with time), discrete variables, and constants. A physical entity has one effort 
variable, which is a continuous variable and the "main" attribute of the physical entity. The continuous 
dynamics of a physical entity is given by a differential equation which is used to compute the value 
of the effort variable. A physical entity can have one or more physical interactions with one or more 
physical entities. The values of the flow variables of these interactions are used in the computation of the 
continuous dynamics of the physical entity. 

A physical interaction (such as the heat flow from the coffee to the room) represents an interaction 
between two physical entities. It consists of one flow variable, a set of attributes, and a continuous dy- 
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Figure 2: Physical system behaviors. 



namics. The flow variable is a continuous variable. The continuous dynamics is described by a differen- 
tial equation. The values of the effort variables from the two physical entities are used in the computation 
of the continuous dynamics of the physical interaction. 

4.1.1 Physical System Behavior 

The basic behavior of physical system components is their continuous behavior. However, there are 
natural phenomena that exhibit a combination of continuous and discrete behaviors. In our method, the 
continuous behavior is performed as long as time advances. When a discrete event must happen, time 
cannot advance until the event has been performed. Figure [2] shows the execution of physical system 
behaviors. 

Approximating continuous behaviors. We use a numerical approach to approximate the continuous 
behaviors of the physical system components by advancing time in discrete time steps, and computing 
the values of the continuous variables at each "visited" point in time. There are two main steps in the 
computation of the continuous behavior which is performed in each time step. The first step is to compute 
the new value of the flow variable of each physical interaction in the system. This computation is based 
on the continuous dynamics of the physical interaction, the values of its attributes, and the "old" values 
(i.e., the values computed at the previous visited point in time) of the effort variables of the two physical 
entities connected to the interaction. In the configuration of physical system components in Figure [T] 
the new values of the flow variables f\-2 to f\ :n of the physical interactions PI\a to PI\ M , respectively, 
are computed. The computation of the new value of f\ : 2 of PI\ : 2, for example, is not only based on the 
attributes of the physical interaction, but also on the "old" values of e\ of PE\ and e2 of PE2. 

The second step is to compute the "new" value of the effort variable of each physical entity in the 
system. Beside the current value of the effort variable, the continuous dynamic of the physical entity, and 
the values of its attributes, the computation also needs the sum of the (new) values of the flow variables 
of all the physical interactions connected to it. For example, in the system in Fig. [T] the new values of the 
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effort variables e\ to e„ of the physical entities PE[ to PE n , respectively, are computed. The computation 
of the new value of e\ of PE\, for example, also depends on the sum of the new values of f\a to f\- n . 

The difference between computing the new values of the effort variables and the flow variables is 
that for the former we need a differential equation solution technique, whereas for the latter we do not 
need such a technique since a flow variable represents a rate. 

Discrete behaviors. The conditions for triggering discrete events depend on the values of the continu- 
ous variables. We consider two kinds of discrete events: one that changes the continuous dynamics of a 
physical entity or interaction (e.g., when the coffee starts boiling, its continuous dynamics is changed), 
and one that changes the value of a discrete variable. The change caused by a discrete event will effect 
the computation of the continuous behavior in the next time step. 

4.2 Applying the Framework for Modeling Thermal Systems 

We can instantiate the above framework to thermal systems as follows: a physical entity is a thermal 
entity, whose effort variable (T) denotes the temperature of the entity and whose continuous dynamics is 
a differential equation defining the heat gained or lost by the entity as its temperature changes. Likewise, 
a physical interaction is a thermal interaction whose flow variable (Q) denotes the heat flow rate. We 
may have different kinds of thermal interactions: conduction, convection, and radiation; their continuous 
dynamics are given by differential equations used for the heat transfer rates. Table [T] shows the thermal 
system components, and Fig. [3] shows a model of our example with a hot cup of coffee in a room. 



Table 1: Thermal system components. 



Components 


Effort/Flow 


Attributes 


Continuous dynamic 


Thermal entity 


Effort: T 


m, c 


f = m 

m-c 


Thermal interaction: conduction 


Flow: Q 


A, k,L 


Q=f-(Tx-T 2 ) 


Thermal interaction: convection 


Flow: Q 


A,h 


Q = h-A-(T l -T 2 ) 


Thermal interaction: radiation 


Flow: Q 


A, e, a 


Q = e-o-A-(T 1 4 -T 2 4 ) 




Figure 3: A cup of coffee in a room with the convection and conduction heat transfers between them. 

5 Modeling Thermal Systems in Real-Time Maude 



This section explains how we can apply the general methods for modeling physical systems given in 
Section|4]to model and execute thermal systems in Real-Time Maude. In particular, Section 5.1 defines 
classes for thermal entities and interactions, Section 5.2 explains how the Euler method (see, e.g., JH) 
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can be used to define in Real-Time Maude an approximation of the continuous dynamics of a thermal 



system, Section 5.3 describes how hybrid behaviors of thermal systems can be modeled, and Section 5.4 
explains how external heat sources can be easily added to our models. 



5.1 Defining Thermal System Components 

In our method, thermal entities (physical entities in thermal systems) can be defined in Real-Time Maude 
as object instances of the following class ThermalEntity: 

class ThermalEntity I heatCap : PosRat, mass : PosRat, 

temp : Rat, tempDisplay : String, 

heat Trans : Rat , mode : CompMode . 

sort CompMode . op default : -> CompMode [ctor] . 

The attribute heatCap denotes the heat capacity and mass the mass of the entity; temp denotes the tem- 
perature, which corresponds to the effort variable in thermal systems; tempDisplay is used to display 
the temperature in a more readable format than as a rational number; finally, heat Trans and mode are 
needed, as we will see later, to implement phase transitions, i.e., to encode the computation mode of a 
thermal entity, which is related to its discrete behavior. 

We define a class for general thermal interactions (physical interactions in thermal systems), and 
three subclasses for the three different heat transfer mechanisms: 

class Thermallnteraction I entityl : Did, entity2 : Qid, area : PosRat, 

Qdot : Rat, QdotDisplay : String . 

class Conduction I thermCond : PosRat, thickness : PosRat . 
class Convection I convCoeff : PosRat . 
class Radiation I emmissiv : PosRat . 

subclass Conduction Convection Radiation < Thermallnteraction . 

In the Thermallnteraction class, Qdot is the attribute corresponding to the heat flow rate Q of the 
thermallnteraction. QdotDisplay is used to display the heat flow rate in a readable format, entityl and 
entity2 are object identifiers of the two thermal entities connected by an object of this class, and area 
is the area of the interaction common for all three interaction types. The remaining interaction-specific 
attributes, e.g., for conductivity the thermal conductivity thermCond of the material and the thickness 
thickness of the material through which the conduction occurs, are specified in the corresponding 
subclasses. 



5.2 Approximating the Continuous Behaviors using the Euler Method 

We use ordinary differential equations for specifying the continuous dynamics of physical entities and 
interactions. For executing the continuous behavior of a physical system by using a discrete-time sam- 
pling strategy, we need an approximation to compute the values of the effort and flow variables. In this 
paper we adapt the Euler method [5], a numerical approach to solve differential equations. Assume that 
the value of the continuous variable y at time t is defined by the differential equation y = f(y,t). Assume 
furthermore that the value of y at time t n is y n . When time advances, the Euler method uses the function 
y n+ i =y n +h- f(y„,t„) to compute the value of y at time point t n+ \ = t„ + h, where h is the time step. 
To implement the Euler method, we define the following functions : 
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• computeQdot that computes the heat flow rate Q of each thermal interaction, 

• computeTemp that computes the temperature T of each thermal entity in the system, 

• sumQdots that computes the sum of the Qs of the thermal interactions of each entity, and 

• ODESol that computes the next value y n+ \ = y n + h*f(y n ,t n ) of a continuous variable y using the 
Euler method. 

We start with the declaration of the following variables: 

vars T Tl T2 QDOT YN FYN QT NEW-Qdot CURR-TEMP NEW-HEAT-TRANS : Rat . 
vars k A L EPSILQN c m h : PosRat . var S : DnOff . 

vars E El E2 TI HG : Oid . vars CONFIG REST : Configuration . 

The function computeQdot computes the heat flow rate Q of each thermal interaction in the system 
according to the laws given in Section [3] and Table [T] The following equation computes the Qdot (given 
by Q = • {Ti — Tj) as explained in Section |3j) of a thermal conduction TI between two thermal entities 
El and E2, and then recursively applies the computeQdot function on the remaining configuration: 

op computeQdot : Configuration ~> Configuration . 
ceq computeQdot ( 

< El : ThermalEntity I temp : Tl > 

< E2 : ThermalEntity I temp : T2 > 

< TI : Conduction I entityl : El, entity2 : E2, thermCond : k, area : A, thickness : L > 
REST) 

< TI : Conduction I Qdot : NEW-Qdot , QdotDisplay : display {NEW-Qdot , precision) > 
computeQdot (< El : ThermalEntity I > < E2 : ThermalEntity I > REST) 

if NEW-Qdot := (k * A * (Tl - T2) ) / L . 

Likewise, the following two equations compute the new value of Qdot of, respectively, thermal 
convections and radiations, and the last equation removes the computeQdot operator when there are no 
more thermal interactions in the remaining configuration: 

ceq computeQdot ( 

< El : ThermalEntity I temp : Tl > 

< E2 : ThermalEntity I temp : T2 > 

< TI : Convection I entityl : El, entity2 : E2, convCoeff : h, area : A > REST) 

< TI : Convection I Qdot : NEW-Qdot , QdotDisplay : display (NEW-Qdot , precision) > 
computeQdot (< El : ThermalEntity I > < E2 : ThermalEntity I > REST) 

if NEW-Qdot := h * A * (Tl - T2) . 

ceq computeQdot ( 

< El : ThermalEntity I temp : Tl > 

< E2 : ThermalEntity I temp : T2 > 

< TI : Radiation I entityl : El, entity2 : E2, emmissiv : EPSILQN, area : A > REST) 

< TI : Radiation I Qdot : NEW-Qdot , QdotDisplay : display (NEW-Qdot , precision) > 
computeQdot (< El : ThermalEntity I > < E2 : ThermalEntity I > REST) 

if NEW-Qdot := EPSILON * stefBolzConst * A * ((Tl ' 4) - (T2 '4)) . 

eq computeQdot (CONFIG) = CONFIG [owise] . 
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The function computeTemp computes the temperature T of each thermal entity in the system, with 
the dynamics given by f = where ^2 is the sum of the Qs of the thermal interactions of the entity 
which is computed by the function sumQdots below. The function uses a numerical method (ODESol) 
which needs the value of the current temperature (T), the time step used, and the continuous dynamics: 

op computeTemp : Configuration ~> Configuration . 

ceq computeTemp (< E : ThermalEntity I heatCap: C, mass : M, temp : T, mode : default > REST) 

<E : ThermalEntity I temp: CURR-TEMP , tempDisplay : display {CURR-TEMP, precision) > 
computeTemp (REST) 
if CURR-TEMP := ODESol (T, timeStep, sumQdots (REST, E) / (M * C))) . 

eq computeTemp (CONFIG) = CONFIG [owise] . 

The function sumQdots computes the sum values of the heat flow rate of all thermal interactions 
connected to a thermal entity. The direction of the heat flow determines whether a thermal entity gains 
or loses the heat. This is represented by multiplying the heat flow rate value by minus one for one part 
of the interaction: 

op sumQdots : Configuration Oid ~> Rat . 

eq sumQdots (< TI : Thermallnteraction I entityl : El, entity2 : E2, Qdot : QDOT > REST, E) 

if (E == El or E == E2) then 

(if E == El then -1 * QDOT + sumQdots (REST, E) else QDOT + sumQdots (REST, E) fi) 
else sumQdots (REST, E) fi . 

eq sumQdots (CONFIG, E) = [owise] . 

The function ODESol uses the Euler method to compute the next value y n+ \ = y n + h* f(y n ,t n ) of a 
continuous variable y, given the current value y n of the variable, the time step h, and the value f(y n ,t n ) 
computed from the continuous dynamics of the variable: 

op ODESol : Rat PosRat Rat -> Rat . 

eq ODESol (YN, H, FYN) = YN + H * FYN . 

As explained in Section[2j a tick rule models the advance of time in a system. We define the following 
tick rule that advances time by one time step and computes the new temperatures for all entities: 

rl [tick] : {CONFIG} => {computeTemp(computeQdot (CONFIG) , timeStep))} in time timeStep . 

5.3 Defining Thermal System with Hybrid Behaviors 

As mentioned in Section |4.1.1| a physical system may also exhibit (at least) two kinds of discrete be- 
haviors. In this paper we focus on discrete events that change the continuous dynamics of an entity or 
interaction. In thermal systems, we have phase transitions between different phases of the entities. For 
the example of water, a discrete event changes the phase from "solid" to "melting," thereby also changing 
the continuous dynamics of the water: the received thermal energy is used for the melting process and 
the temperature of the water does not change during melting. The following changes need to be made to 
our model to accommodate such hybrid behaviors in physical systems: 

1. The basic ThermalEntity class must be extended to keep track of the phases of the entity. 

2. Instantaneous rewrite rules modeling the discrete change from one phase to another must be added. 
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3. It must be ensured that the above rules are applied at the right time. 

4. The functions computing the values of the continuous variables must take into account the different 
continuous dynamics in different phases. 

Extending the basic class. Let's say we want to define a thermal entity representing a water substance 
with three phases: solid, liquid, and gas. Between these main phases, the water is either melting, freezing, 
evaporating, or condensing. Figure |4] shows the various phases (and their corresponding continuous 
dynamics) and the transitions between them. 




Figure 4: The phases of the water entity. 



The class WaterEntity extending the base class ThermalEntity can be defined as follows: 

class WaterEntity I phase : MatterState, heat Trans : Rat, heatTransDisplay : String . 
subclass WaterEntity < ThermalEntity . 

sort MatterState . 

ops liquid solid gas melting evaporating condensing freezing : -> MatterState . 
op phaseChange : -> CompMode . 

The new attribute phase denotes the current phase of the water substance, and heat Trans denotes the 
latent heat of the water in the transitions between the three main phases. Notice that the continuous 
dynamics of the temperature is the same in all the three main phases of water, whereas the temperature 
does not change in-between these phases. In addition to the def ault computation mode defined above, 
we also add a new computation mode phaseChange to denote these intermediate phases of water, so 
that the computeTemp function is defined correctly. 



Rewrite rules denning discrete change. A phase change is modeled by an instantaneous rewrite rule. 
We show two of the eight such rules (corresponding to the arrows in Fig. [4]) for water: 

crl [solid-to-melting] : 

< E : WaterEntity I temp : T, phase : solid > 
=> 

< E : WaterEntity I phase : melting, mode : phaseChange, heatTrans : > 
if T >= . 



crl [melting-to-liquid] : 

< E : WaterEntity I phase : melting, heatTrans : QT, mass : m > 
=> 

< E : WaterEntity I phase : liquid, mode : default > 
if QT / m >= latentHeatFus . 
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The change from a "main" phase to a "transitional" phase happens when the temperature reaches a given 
value, whereas a transition the other way happens when the value of the heat accumulated during the 
transitional phase divided by the mass of the entity reaches some value (such as the latent heat fusion). 

We ensure that the above rules are applied when they should by disabling the tick rule when one of 
the above rules is enabled. This is done by defining a function timeCanAdvance on the state, so that 
timeCanAdvance holds only if no such instantaneous rule needs to be applied in the current state. We 
then modify the tick rule so that it can only be applied to states where timeCanAdvance holds. Again, 
we only show two of the equations defining timeCanAdvance and the modified tick rule: 

eq timeCanAdvance (< E : WaterEntity I phase : solid, temp : T >) = T < . 
eq timeCanAdvance (< E : WaterEntity I phase : melting, heatTrans : QT, mass : m >) = 
QT / m < latentHeatFus . 

crl [tick] : {CONFIG} => {computeTemp (computeQdot (CONFIG) , timeStep)} in time timeStep 
if timeCanAdvance (CONFIG) . 



The change of continuous dynamics caused by a phase change means that there is more than one 
possible continuous dynamics for a thermal system component. We need to modify the function used 
to compute the corresponding continuous variable. For the water entity, we need to modify the function 
computeTemp, so that when the entity is in mode phaseChange (as opposed to the default mode, 
where the equation in Section 5.2 applies), then the temperature does not change, but the accumulated 
heat is increased/decreased: 



ceq computeTemp (< E : WaterEntity I mode : phaseTrans , heatTrans : QT > REST) = 
< E : WaterEntity I heatTrans : NEW-HEAT-TRANS, 

heatTransDisplay : display (NEW-HEA T-TRANS , precision) > 

computeTemp (REST) 
if NEW-HEAT-TRANS := ODESo I ( QT, timeStep, sumQdots (REST, E)) . 



5.4 External Heat Sources 

It is sometimes convenient to be able to inject heat into a thermal entity. For example, we may want to 
add a boiler to our coffee system that adds heat to the cup of coffee. In our method, such heat generators 
can be defined in Real-Time Maude as object instances of the following class HeatGenerator: 

class HeatGenerator I Qdot : Rat, entity : Oid . 

The attribute Qdot denotes the heat flow rate provided by the heat generator; entity is the object 
identifier of the thermal entity connected to the heat generator. 

Adding a heat generator to a thermal entity means that beside the heat flows from the interactions with 
other thermal entities, the heat flow from the heat generator influences the change of the temperature of 
the entity. We need to modify the function used to compute the sum of heat flow rate values of a thermal 
entity to include the heat flow rate from the connected heat generator: 

eq SumQdots((< HG : HeatGenerator I Qdot : QDOT, entity : El > REST), E) = 
if (E == El) then QDOT + SumQdots (REST, E) else SumQdots (REST, E) fi . 
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Table 2: Initial values for object constant parameters (n = y). 



room 



cup of coffee 



airDensity 

roomVolume 

roomMass 

roomHC 

h 



air density Fj 

volume 

mass 

heat capacity 



12 
10 



64 m 



384 
5 

105 
10 
20 
1000 



kg/nr' 
3 



cupRadius 

cupHeight 

cupCircum 

cupBaseArea 

cupSideArea 

cupThickness 

waterDensity 

cof f eeVolume 

cof f eeMass 

coffeeHC 

k 



inner radius 
height 

circumference 
circle area 
side area 
cup thickness 
water density 
volume 
mass 

heat capacity 



3. 



4 

100 
9 

100 
44 
175 

22 
4375 

99 
4375 
1 

200 
1000 



218750 
396 
875 
42 
10 

15 
10000 



Free convection of gases is between 2-25 
c Thermal conductivity of porcelain 



6 Case Studies 



m 

kg/nr 1 

2 

m 
kg 

kJ/kg 

kW 
m-°C 



fl At sea level and 20°C, air has a density of approximately 1. 2kg (http://en.wikipedia.org/wiki/Density_of_air). 



This section shows how our method can be used to model and analyze some simple thermal systems. 



In particular, in Section 6.1 we model and analyze a cup of hot coffee in a room, with two kinds of 
interactions {conduction and convection). In Section [6T2] we add a heater that sends constant heat to the 



cup of coffee. Finally, in Section 6.3 we add a "smart" heater that switches itself on or off in order to 
keep the temperature of the coffee in a desired interval. It is worth noticing that, using the definitions 
in the previous section, the definition of the first two models reduces to defining appropriate initial 
states. Table [2] shows the constants used for modeling the room and the cup of coffee. The analyses 
are performed on an AMD Athlon(tm) 64 Processor 3500+ with 2GB of RAM, using a timeStep of 1, 
unless something else is specified. 



6.1 Case Study 1: A Cup of Coffee in a Room 

In our first case study we model the simple coffee system in Fig. [3] in Section |4~2| by defining the corre- 
sponding initial state as follows: 

op csl : -> GlobalSystem . 
eq csl = 

{< coffee : ThermalEntity I heatCap : coffeeHC, mass : cof f eeMass, heatFlow : 0, 

temp : 70, tempDisplay : "", mode : default, heatTrans : > 

< room : ThermalEntity | heatCap : roomHC, mass : roomMass, temp : 20, mode : default, 

tempDisplay : "", heatFlow : 0, heatTrans : > 

< crConduct : Conduction I entityl : coffee, entity2 : room, thermCond : k, Qdot : 0, 

area : condArea, thickness : cupThickness, QdotDisplay : "" > 

< crConvect : Convection I entityl : coffee, entity2 : room, convCoeff : h, 

area : convArea, Qdot : 0, QdotDisplay : "" >} . 
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where coffeeHC, coffeeMass, ... are constants of sort Rat with values as shown in Table|2] The 
constant condArea is the conduction area, computed as the sum of the coffee circle area on the base of 
the cup and the coffee cylindrical area on the lateral surface of the cup. The constant convArea is the 
convection area, equal to the coffee circle area, representing the top surface of the coffee. 

We can simulate the behavior from this initial state up to a certain time, using Real-Time Maude's 
timed rewrite command: 

Maude> (trew csl in time <= 1000 .) 



result ClockedSystem : 

{< crConduct : Conduction 

< crConvect : Convection 

< coffee : ThermalEntity 

< room : ThermalEntity I 



I QdotDisplay 
I QdotDisplay 
I TempDisplay 
TempDisplay : 1 



: "0.0044820616", 
: "0.0000543280", 
: "21.6767974687" 
21.1390469168", 



>> in time 1000 



We show only the display attributes of the objects. The simulation takes about a minute and a half. 
The result is as expected: in 1000 seconds, the temperature of the coffee is sensibly decreased, reaching 
almost that of the room, whose temperature has increased slightly due to the heat from the coffee. 

Moreover, we can collect the values for coffee and room temperatures at each time step during the 
simulation. This has been done by adding a "collector" object to the initial configuration, which, at each 
tick rule application, records the current time and temperature values. This allows us to plot the coffee 



and room temperatures as a function over time, as shown in Figure i '&) 




Time Time 



(a) Case study 1 : simulation (b) Case study 2: simulation 

Figure 5: Simulation plots for the case studies 1 and 2. 



Suppose we would like to check how long it takes for the coffee and room to reach about the same 
temperature, with a given tolerance. This can be done by searching for one reachable state where the 
difference of the two temperatures is less than the given tolerance. This kind of search can be performed 
in Real-Time Maude without any restriction on time, since we can always reach such a desired state 
within a finite amount of time: 

Maude> (tsearch [1] csl =>* 

{< coffee : ThermalEntity I temp : Tcoffee:Rat > 
< room : ThermalEntity / temp : Troom:Rat > C: Configuration} 

such that (abs (Tcoffee:Rat - TroomiRat) <= 1/1000 ) with no time limit .) 
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Figure 6: Case study 2: Modeling a cup of coffee in a room with a manual coffee heater. 

It takes about 17 minutes for Real-Time Maude to find a solution, which indicates that it takes around 
2340 seconds for the room and coffee to reach the same temperature. 

6.2 Case Study 2: Adding a Boiler to the Cup of Coffee 

We extend our first case study by adding a heater that is always on and therefore sends a constant flow of 
heat to the cup of coffee, as illustrated in Fig. [6] Furthermore, we start with an initial coffee temperature 
of — 10°C, so that this example shows how the coffee entity exhibits a hybrid behavior, when given a 
sufficient constant heat, going through all three main phases: solid, liquid, and gas. 

We now declare the coffee object to be a WaterEntity object, and add an object boiler of class 



HeatGenerator. The other objects (room, crConduct, and crConvect) are defined as in Section 6.1 



op cs2 : -> GlobalSystem . 
eq cs2 = 

{< coffee : WaterEntity I heatCap : coffeeHC, mass : coffeeMass, temp : -10, 

tempDisplay : "", heatFlow : 0, mode : default, phase : solid, 
heat Trans : 0, heatTransDisplay : "" > 

< boiler : HeatGeneration I entity : coffee, Qdot : 15/10 >} . 

As for the first case study, we can simulate the above system and plot the room and coffee tempera- 



tures as functions of time. Figure f b) shows the simulation up to time 1500. We notice that the coffee 
temperature remains constant during the transition phases (melting and evaporation), while it increases 
when the coffee is in liquid, solid, or gaseous state. The fact that the melting temperature is somewhat 
greater than °C is a consequence of the approximation in the computation of the continuous behavior 
and the chosen timeStep: a bigger timeStep will lead to a coarser numerical approximation in the 
Euler algorithm. The room tends to become warmer, compared to the first case study; this is due to the 
presence of the heater, which, through the coffee, indirectly transfers heat to the room. 

If we would like to know how long it takes for the iced coffee to start to melt, we can use Real-Time 
Maude's find earliest command to find the first reachable state such that the coffee is melting: 

Maude > (find earliest cs2 =>* {C: Configuration < coffee : WaterEntity I phase : melting >} .) 
In 224ms Real-Time Maude returns a solution showing that the iced coffee starts to melt after 22 seconds. 

6.3 Case Study 3: Keeping the Coffee Warm 

We keep the set-up of the second case study, but now we define a more sophisticated heater, namely, 
one that senses the coffee temperature at each time step and tries to keep the temperature of the coffee 
between 70°C and 80°C by turning itself on and off. We need to modify the previous model by: 
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• adding two instantaneous rules: one for turning the heater on when the coffee temperature is below 
70°C, and one for turning the heater off when the coffee is warmer than 80°C, and 

• adding a condition to the tick rule so that time does not advance when one of the above rules should 
be taken. 

For convenience, we define a subclass SmartHeater of the class HeatGenerator for such smart heaters: 

class SmartHeater I status : OnOff, lowTemp : Rat, highTemp : Rat . 
subclass SmartHeater < HeatGenerator . 



sort OnOff . ops on off : -> OnOff . 

The rules for turning the heater on and off are immediate: 

crl [turnOff] : 

< E : ThermalEntity I temp : Tl > 

< HG : SmartHeater I entity : E, status : on, highTemp : T2 > 
=> 

< E : ThermalEntity I > 

< HG : SmartHeater I status : off , Qdot : > 
if Tl >= T2 . 



crl [turnOn] : 

< E : ThermalEntity I temp : Tl > 

< HG : SmartHeater I entity : E, status : off , lowTemp : T2 > 
=> 

< E : ThermalEntity I > 

< HG : SmartHeater I status : on, Qdot : 15/10 > 
if Tl <= T2 . 

We also define a function heatersOK, which holds in a state when none of the above two rules can 
be applied, and modify the tick rule to take this into account: 

op heatersOK : Configuration -> Bool [frozen (1)] . 

eq heatersOK (< HG : SmartHeater I entity : E, status : S, lowTemp : Tl, highTemp : T2 > 
< E : ThermalEntity I temp : T > REST) 
= ((S == on and T < T2) or (S == off and T > Tl)) 
and heatersOK ( < E : ThermalEntity I > REST) . 
eq heatersOK (CONFIG) = true [owise] . 

crl [tick] : {CONFIG} => {computeTemp (computeQdot (CONFIG) , timeStep)} in time timeStep 
if timeCanAdvance (CONFIG) and heatersOK (CONFIG) . 

The initial state is as in the second case study, with the exception that the "dumb" heater has been 
replaced by a smart heater: 

op cs3 : -> GlobalSystem . 
eq cs3 = 

{< coffee : WaterEntity I temp : -20, phase : liquid, ... > 
< coffeeHeater : SmartHeater I entity : coffee, status : off, Qdot : 0, lowTemp : 70, 

highTemp : 80 > 

...} . 
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Figure 7: Case study 3: Simulation. 

Figure [7] shows the coffee and room temperatures as functions of the elapsed time in a Real-Time 
Maude simulation of this system. 

Finally, we use time -bounded LTL model checking to analyze the stability property that once the 
temperature of the coffee has reached 69.5°C, its temperature will remain in the interval between 69.5°C 
and 80.5°C (we set the interval as 69.5 to 80.5 instead of 70 to 80 for taking into account the imprecision 
from the approximation using the Euler method). We define the atomic proposition temp-ok to hold in 
all states where the coffee temperature is between 69.5°C and 80.5°C: 

op temp-ok : -> Prop [ctor] . 

ceq {REST < coffee : WaterEntity I temp : T >} |= temp-ok = (T >= 139/2 and T <= 161/2) . 

The property can be checked by executing using time-bounded model checking: 
Maude> (me cs3 l=t [] (temp-ok -> [] temp-ok) in time <= 1500 .) 

Result Bool : 
true 

7 Related Work 

Modelica [3] is an object-oriented language for modeling physical systems where hybrid differential 
algebraic equations model the continuous dynamics. The language supports the specification of linear 
and non-linear equations. Since the language does not have a formal semantics, the precise meaning of 
a model depends on the simulation tool used. Furthermore, there are no reachability and temporal logic 
analysis tools for Modelica. 

Hy Visual (8) is an actor-based tool environment for the simulation of continuous and hybrid systems. 
Hybrid behaviors are specified as finite state machines in which a state can be refined into a dynamic sys- 
tem represented as ordinary differential equations. Hy Visual does not provide model checking analyses. 

OO-DPT [16l is an object-oriented approach for modeling hybrid supervisory systems. This ap- 
proach proposes a Petri-net-based formal technique named object-oriented differential predicate transi- 
tion net (OO-DPT net). This technique defines an interface between differential equation systems and 
Petri nets to model hybrid behaviors. Reachability and safety properties analysis is the target of this 
approach, but as of yet, there is no tool support for OO-DPT. 



coffee 

room 
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HyTech Hi is a tool for modeling and verifying hybrid systems represented as linear hybrid automata. 
The discrete behaviors are represented as a set of locations and discrete transitions between them where 
the continuous dynamics is defined in each location. The continuous dynamics is specified as polyhedral 
differential inclusions. The restrictions related to the formalism used by the tool confines the representa- 
tion of continuous dynamics to the ones usually used in modeling physical systems. The main difference 
between HyTech and Real-Time Maude is the expressiveness and generality of the formalism, where in 
Real-Time Maude all kinds of data types, functions, communication models, dynamic object creation, 
etc., can be specified. 

SHYMaude (Stochastic Hybrid Maude) |[T0ll is a modeling language for object-based stochastic 
hybrid systems that extends the PMaude [7 ] tool for probabilistic rewrite theories. Stochastic hybrid 
systems consist of distributed stochastic hybrid objects that interact with each other by asynchronous 
message passing. The continuous behaviors are governed by stochastic differential equations, and the 
discrete changes are probabilistic. For simulation, the continuous dynamics of stochastic differential 
equations is approximated by using the Euler-Maruyama method, a generalization of Euler method to 
approximate numerical solution of stochastic differential equations. For formal analyses, the interface of 
the statistical model checking tool VeStA [ 15] to Maude system is used. Despite the obvious similarities 
that both approaches are object-oriented approaches to hybrid systems in extensions of Maude, our work 
and the SHYMaude work are in fact quite different. In [10] the hybrid objects are autonomous and only 
interact with other hybrid objects by message passing, whereas in our work, the main focus is on the 
physical interaction (such as heat flow) between the hybrid objects. 

8 Concluding Remarks 

We have presented a general object-based method for formally modeling and analyzing thermal systems 
in Real-Time Maude. In contrast to most other approaches, we also focus on the physical interactions 
between different physical components. We explain how the Euler method can be adapted to approximate 
continuous behaviors, and have illustrated our method on three variations of a simple coffee-and-room 
system with realistic parameters. 

This work should be seen as a first exploration into how hybrid systems can be modeled and analyzed 
in Real-Time Maude. Although we believe that Real-Time Maude - because of its expressiveness and 
generality, support for objects, user-definable data types, different communication models, etc. - should 
be a suitable candidate for modeling and analyzing advanced hybrid systems; this of course has to be 
validated by applying our techniques on advanced hybrid systems. 

Euler is a simple numerical algorithm for solving ordinary differential equations (ODE). It could be 
interesting to implement in Real-Time Maude some other numerical approaches for solving ODE (e.g. 
the midpoint method, the linear multi-step method, or the Runge-Kutta method) that offer more precise 
results than Euler, probably at the cost of computational efficiency. Another point to explore is the 
possibility of using the other time sampling strategies offered by Real-Time Maude, like the maximum 
time advance strategy; applying this strategy to continuous behavior would require an execution strategy 
that can "predict" the maximum time that the system can advance, before some instantaneous rule has to 
be taken, possibly using an advanced numerical method for ODE with dynamic step. 
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